SUBROUTINE VDFDIFH1S(&
 & KIDIA  , KFDIA  , KLON   , KLEV   , KTOP   , KTILES , KTVL   , KTVH , KLEVSN,&
 & PSSDP2 , PTMST  , PFSH1D , PFLH1D , LDFLUX1D, &
 & PFRTI  , PCIL, PSSRFLTI,PSLRFL , PEMIS  , PEVAPSNW, &
 & PHLICE , PTLICE , PTLWML , &
 & PSLM1  , PTM1   , PQM1   , PQTM1  , PAPHM1 , &
 & PCFH   , PCFHTI , PCFQTI , PMFLX  , PSLUH  , PQTUH  , &
 & PTDIF  , PQDIF  , PCPTSTI, PQSTI  , PCAIRTI, PCSATTI, &
 & PCPTSTIU,PCAIRTIU,PCSATTIU,PTSRF,PLAMSK,PSNM,PRSN, &
 & PDQSTI , PTSKTI , PTSKRAD, PTSM1M , &
 & PTSNOW , PTICE  , PSST, &
 & PTSKTIP1, PSLGE  , PTE    , PQTE, &
 & PJQ,PSSH,PSLH,PSTR,PG0,PJQU)  
! (C) Copyright 1989- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!     ------------------------------------------------------------------

!**   *VDFDIFH* - DOES THE IMPLICIT CALCULATION FOR DIFFUSION OF S. L.

!     DERIVED FROM VDIFF (CY34) BY
!     A.C.M. BELJAARS       E.C.M.W.F.    10-11-89

!     OBUKHOV-L UPDATE      ACMB          26/03/90.
!     SKIN T CLEANING       P. VITERBO    15-11-96.
!     TILE BOUNDARY COND.   ACMB          20-11-98.
!     Surface DDH for TILES P. Viterbo    17-05-2000.
!     New tile coupling, 
!     DDH moved to VDFMAIN  A. Beljaars   2-05-2003.
!     Mass flux terms,
!     Flux b.c. for SCM,
!     Moist generalization  A. Beljaars/M. Ko"hler 3-12-2004. 
!     Removed option for linearized    P. Lopez  02/06/2005
!     physics (now called separately)   

!     PURPOSE
!     -------

!     SOLVE TRIDIAGONAL MATRICES FOR DIFFUSION OF DRY STATIC ENERGY
!     AND MOISTURE; IN SO DOING, IT ALSO SOLVES THE SKIN TEMPERATURE
!     EQUATION.

!     INTERFACE
!     ---------

!     *VDFDIFH* IS CALLED BY *VDFMAIN*

!     INPUT PARAMETERS (INTEGER):

!     *KIDIA*        START POINT
!     *KFDIA*        END POINT
!     *KLEV*         NUMBER OF LEVELS
!     *KLON*         NUMBER OF GRID POINTS PER PACKET
!     *KTOP*         INDEX FOR BOUNDARY LAYER TOP
!     *KTVL*         DOMINANT LOW VEGETATION TYPE
!     *KTVH*         DOMINANT HIGH VEGETATION TYPE

!     INPUT PARAMETERS (REAL):

!     *PTMST*        DOUBLE TIME STEP (SINGLE AT 1TH STEP)
!     *PFRTI*        FRACTION OF SURFACE AREA COVERED BY TILES
!     *PSLM1*        GENERALIZED LIQUID WATER STATIC ENERGY    AT T-1
!                    (NOTE: In lin/adj physics = Dry static energy)
!     *PTM1*         TEMPERATURE AT T-1
!     *PQM1*         SPECIFIC HUMIDITY      AT T-1
!     *PQTM1*        SPECIFIC TOTAL WATER   AT T-1
!     *PAPHM1*       PRESSURE AT T-1
!     *PCFH*         PROP. TO EXCH. COEFF. FOR HEAT (C,K-STAR IN DOC.)
!     *PCFHTI*       IDEM FOR HEAT (SURFACE LAYER ONLY)
!     *PCFQTI*       IDEM FOR MOISTURE (SURFACE LAYER ONLY)
!     *PMFLX*        MASSFLUX
!     *PSLUH*        UPDRAFT GENERALIZED LIQUID WATER STATIC ENERGY AT HALF LEVEL
!     *PQTUH*        UPDRAFT SPECIFIC TOTAL WATER AT HALF LEVEL
!     *PCPTSTI*      DRY STATIC ENRGY AT SURFACE
!     *PQSTI*        SATURATION Q AT SURFACE
!     *PCAIRTI*      MULTIPLICATION FACTOR FOR Q AT LOWEST MODEL LEVEL
!                    FOR SURFACE FLUX COMPUTATION
!     *PCSATTI*      MULTIPLICATION FACTOR FOR QS AT SURFACE
!                    FOR SURFACE FLUX COMPUTATION
!     *PCPTSTIU*     AS PCPTSTI FOR UNSTRESSED EVAPORARTION FROM LOW VEGET.
!     *PCAIRTIU*     AS PCAIRTI FOR UNSTRESSED EVAPORARTION FROM LOW VEGET.
!     *PCSATTIU*     AS PCSATTI FOR UNSTRESSED EVAPORARTION FROM LOW VEGET.
!     *PDQSTI*       D/DT (PQS)
!     *PSSRFLTI*     NET SOLAR RADIATION AT THE SURFACE, FOR EACH TILE
!     *PSLRFL*       NET THERMAL RADIATION AT THE SURFACE
!     *PEMIS*        MODEL SURFACE LONGWAVE EMISSIVITY
!     *PEVAPSNW*     EVAPORATION FROM SNOW UNDER FOREST
!     *PTSKTI*       SKIN TEMPERATURE AT T-1
!     *PTSKRAD*      SKIN TEMPERATURE OF LATEST FULL RADIATION TIMESTEP
!     *PTSM1M*       TOP SOIL LAYER TEMPERATURE
!     *PTSNOW*       SNOW TEMPERATURE   
!     *PTICE*        ICE TEMPERATURE (TOP SLAB)
!     *PSST*         (OPEN) SEA SURFACE TEMPERATURE
!     *PSLGE*        GENERALIZED DRY STATIC ENERGY TENDENCY
!                    (NOTE: In lin/adj physics = Temperature tendency)
!     *PQTE*         TOTAL WATER TENDENCY
!                    (NOTE: In lin/adj physics = Humidity tendency)
!     *ZTSRF*        Surface temperature under each tile
!     *ZLAMSK*       Skin layer conductivity for each tile 

!     OUTPUT PARAMETERS (REAL):

!     *PTDIF*        SLG-DOUBLE-TILDE DEVIDED BY ALFA
!     *PQDIF*        QT-DOUBLE-TILDE DEVIDED BY ALFA
!     *PTSKTIP1*     SKIN TEMPERATURE AT T+1
!     *PJQ*          Surface moisture flux                      (kg/m2s)
!     *PSSH*         Surface sensible heat flux                 (W/m2)
!     *PSLH*         Surface latent heat flux                   (W/m2)
!     *PSTR*         Surface net thermal radiation              (W/m2)
!     *PG0*          Surface ground heat flux (solar radiation  (W/m2)
!                    leakage is not included in this term)
!     *PJQU*         Surface moisture flux unstressed low veg   (kg/m2s)

!     Additional parameters for flux boundary condtion (in 1D model):

!     *LDFLUX1D*     If .TRUE. flux boundary condtion is used 
!     *PFSH1D*       Specified sensible heat flux (W/m2)
!     *PFLH1D*       Specified latent heat flux (W/m2)

!     METHOD
!     ------

!     *LU*-DECOMPOSITION (DOWNWARD SCAN), FOLLOWED BY SKIN-TEMPERATURE
!     SOLVER, AND BACK SUBSTITUTION (UPWARD SCAN).

!     EXTERNALS.
!     ----------

!     *VDFDIFH* CALLS:
!         *SURFSEB*

!     ------------------------------------------------------------------

USE PARKIND1  ,ONLY : JPIM     ,JPRB       ,JPRD
USE YOMHOOK   ,ONLY : LHOOK    ,DR_HOOK, JPHOOK
USE YOEVDF   , ONLY : RVDIFTS
USE YOMCST   , ONLY : RCPD     ,RG     ,RLSTT  ,RLVTT,RTT
USE YOETHF   , ONLY : RVTMP2
USE YOEPHY   , ONLY : LEFLAKE , LESNML, LEURBAN
USE YOMDPHY  , ONLY : YDSURF
IMPLICIT NONE

#include "surfseb.h"


!*         0.1    GLOBAL VARIABLES

INTEGER(KIND=JPIM),INTENT(IN)    :: KLON 
INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTILES 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTVL(KLON) 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTVH(KLON) 
INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTOP
INTEGER(KIND=JPIM),INTENT(IN)    :: KLEVSN
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSSDP2(:,:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMST 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PFSH1D(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PFLH1D(KLON) 
LOGICAL           ,INTENT(IN)    :: LDFLUX1D 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PFRTI(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCIL(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSSRFLTI(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSLRFL(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PEMIS(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PEVAPSNW(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PHLICE(KLON)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTLICE(KLON)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTLWML(KLON)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSLM1(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTM1(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQM1(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQTM1(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPHM1(KLON,0:KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCFH(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCFHTI(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCFQTI(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PMFLX(KLON,0:KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSLUH(KLON,0:KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQTUH(KLON,0:KLEV) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTDIF(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PQDIF(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCPTSTI(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQSTI(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCAIRTI(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCSATTI(KLON,KTILES)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCPTSTIU(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCAIRTIU(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCSATTIU(KLON,KTILES)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSRF(KLON,KTILES)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PLAMSK(KLON,KTILES)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PDQSTI(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSKTI(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSKRAD(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSM1M(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSNOW(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTICE(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSST(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSNM(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRSN(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTSKTIP1(KLON,KTILES)  
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSLGE(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTE(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQTE(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PJQ(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSSH(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSLH(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSTR(KLON,KTILES) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PG0(KLON,KTILES)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PJQU(KLON,KTILES)

!*         0.2    LOCAL VARIABLES

REAL(KIND=JPRB) ::    ZAA(KLON,KLEV) ,ZBB(KLON,KLEV) ,ZCC(KLON,KLEV) ,&
                    & ZTYY(KLON,KLEV),ZQYY(KLON,KLEV),ZGAM(KLON,KLEV),&
                    & Z1DP(KLON,KLEV)
REAL(KIND=JPRB) ::    Z1BET(KLON)    ,&
                    & ZAQL(KLON)     ,ZBQL(KLON)     ,ZASL(KLON)     ,&
                    & ZBSL(KLON)     ,ZSL(KLON)      ,ZQL(KLON)
REAL(KIND=JPRB) ::    ZTSRF(KLON,KTILES)  ,ZTSRFU(KLON,KTILES),ZRHOCHU(KLON,KTILES)  ,& 
                    & ZRHOCQU(KLON,KTILES),ZJS(KLON,KTILES)      ,&
                    & ZSSK(KLON,KTILES)   ,ZTSK(KLON,KTILES) , ZLAMSK(KLON,KTILES)

INTEGER(KIND=JPIM) :: JK, JL, JT

REAL(KIND=JPRB) ::    ZQDP, ZQSP1, ZTPFAC2, ZTPFAC3,&
                    & ZCSNQ, ZCSNS, ZCOEF1,ZCOEF2
REAL(KIND=JPHOOK) ::    ZHOOK_HANDLE

REAL(KIND=JPRB) ::    ZEPSILON

REAL(KIND=JPRB) ::    ZTHKICE(KLON),ZSNTICE(KLON)
LOGICAL         ::    LNEMOLIMTHK,LSICOUP

REAL(KIND=JPRB) :: ZPCAIRTIU(KLON,KTILES),ZPCSATTIU(KLON,KTILES),ZPCPTSTIU(KLON,KTILES),&
                   ZPJQU(KLON,KTILES)
!     ------------------------------------------------------------------

IF (LHOOK) CALL DR_HOOK('VDFDIFH1S',0,ZHOOK_HANDLE)

ZTPFAC2=1.0_JPRB/RVDIFTS
ZTPFAC3=1-ZTPFAC2
ZEPSILON=10._JPRB*EPSILON(ZEPSILON)


!1s !     ------------------------------------------------------------------
!1s 
!1s !*         1.     FULL MODEL PHYSICS WITH MOIST MASS FLUX PBL
!1s !                 -------------------------------------------
!1s 
!1s !*         1.0    SETTING OF THE MATRIX A, B AND C.
!1s 
!1s DO JK=KTOP+1,KLEV-1
!1s   DO JL=KIDIA,KFDIA
!1s     Z1DP(JL,JK)=1.0_JPRB/(PAPHM1(JL,JK)-PAPHM1(JL,JK-1))
!1s     ZAA(JL,JK) =(-PCFH(JL,JK-1)-0.5_JPRB*PMFLX(JL,JK-1))*Z1DP(JL,JK)
!1s     ZCC(JL,JK) =(-PCFH(JL,JK)  +0.5_JPRB*PMFLX(JL,JK)  )*Z1DP(JL,JK)
!1s     ZBB(JL,JK) =1.0_JPRB+(PCFH(JL,JK-1)+PCFH(JL,JK)&
!1s      & -0.5_JPRB*(PMFLX(JL,JK-1)-PMFLX(JL,JK)))*Z1DP(JL,JK)  
!1s   ENDDO
!1s ENDDO
!1s 
!1s !          1.0a   THE SURFACE BOUNDARY CONDITION
!1s 
!1s DO JL=KIDIA,KFDIA
!1s   Z1DP(JL,KLEV)=1.0_JPRB/(PAPHM1(JL,KLEV)-PAPHM1(JL,KLEV-1))
!1s   ZCC(JL,KLEV) =0.0_JPRB
!1s   ZAA(JL,KLEV) =        (-PCFH(JL,KLEV-1)-0.5_JPRB*PMFLX(JL,KLEV-1))*Z1DP(JL,KLEV)
!1s   ZBB(JL,KLEV) =1.0_JPRB+(PCFH(JL,KLEV-1)-0.5_JPRB*PMFLX(JL,KLEV-1))*Z1DP(JL,KLEV)  
!1s ENDDO
!1s 
!1s !          1.0b   THE TOP BOUNDARY CONDITION    
!1s 
!1s   DO JL=KIDIA,KFDIA
!1s     Z1DP(JL,KTOP)=1.0_JPRB/(PAPHM1(JL,KTOP)-PAPHM1(JL,KTOP-1))
!1s     ZAA(JL,KTOP) =0.0_JPRB
!1s     ZCC(JL,KTOP) =         (-PCFH(JL,KTOP)+0.5_JPRB*PMFLX(JL,KTOP))*Z1DP(JL,KTOP)
!1s     ZBB(JL,KTOP) =1.0_JPRB+( PCFH(JL,KTOP)+0.5_JPRB*PMFLX(JL,KTOP))*Z1DP(JL,KTOP)
!1s   ENDDO
!1s 
!1s !*         1.1    SETTING OF RIGHT HAND SIDES.
!1s 
!1s DO JK=KTOP+1,KLEV-1
!1s   DO JL=KIDIA,KFDIA
!1s     ZTYY(JL,JK) = ZTPFAC2 * PSLM1(JL,JK) &
!1s      & + PTMST * PSLGE(JL,JK) &
!1s      & + ZTPFAC2 *(PMFLX(JL,JK)*PSLUH(JL,JK)-PMFLX(JL,JK-1)*PSLUH(JL,JK-1)) *Z1DP(JL,JK)  
!1s     ZQYY(JL,JK) = ZTPFAC2 * PQTM1(JL,JK) &
!1s      & + PTMST * PQTE(JL,JK) &
!1s      & + ZTPFAC2 *(PMFLX(JL,JK)*PQTUH(JL,JK)-PMFLX(JL,JK-1)*PQTUH(JL,JK-1)) *Z1DP(JL,JK)  
!1s   ENDDO
!1s ENDDO
!1s 
!1s !          1.1a   SURFACE
!1s 
!1s JK=KLEV
!1s DO JL=KIDIA,KFDIA
!1s   ZTYY(JL,JK) = ZTPFAC2 * PSLM1(JL,JK) &
!1s    & + PTMST * PSLGE(JL,JK) &
!1s    & + ZTPFAC2 *(-PMFLX(JL,JK-1)*PSLUH(JL,JK-1)) *Z1DP(JL,JK)  
!1s   ZQYY(JL,JK) = ZTPFAC2 * PQTM1(JL,JK) &
!1s    & + PTMST * PQTE(JL,JK) &
!1s    & + ZTPFAC2 *(-PMFLX(JL,JK-1)*PQTUH(JL,JK-1)) *Z1DP(JL,JK)  
!1s ENDDO
!1s 
!1s !          1.1b   TOP
!1s 
!1s JK=KTOP
!1s DO JL=KIDIA,KFDIA
!1s   ZTYY(JL,JK) = ZTPFAC2 * PSLM1(JL,JK) &
!1s    & + PTMST * PSLGE(JL,JK) &
!1s    & + ZTPFAC2 *(PMFLX(JL,JK)*PSLUH(JL,JK)) *Z1DP(JL,JK)  
!1s   ZQYY(JL,JK) = ZTPFAC2 * PQTM1(JL,JK) &
!1s    & + PTMST * PQTE(JL,JK) &
!1s    & + ZTPFAC2 *(PMFLX(JL,JK)*PQTUH(JL,JK)) *Z1DP(JL,JK)  
!1s ENDDO
!1s 
!1s !*         1.2    ADD MOISTURE FLUX FROM SNOW FROM TILE 7 AS EXPLICIT TERM
!1s 
!1s JK=KLEV
!1s DO JL=KIDIA,KFDIA
!1s   ZCSNQ=RG*PTMST*PFRTI(JL,7)*PEVAPSNW(JL)*Z1DP(JL,JK)
!1s   ZCSNS=RCPD*RVTMP2*PTSKTI(JL,7)*ZCSNQ
!1s 
!1s   ZTYY(JL,JK)=ZTYY(JL,JK)-ZCSNS
!1s   ZQYY(JL,JK)=ZQYY(JL,JK)-ZCSNQ
!1s ENDDO
!1s 
!1s !*         1.3    TOP LAYER ELIMINATION.
!1s 
!1s DO JL=KIDIA,KFDIA
!1s   Z1BET(JL)=1.0_JPRB/ZBB(JL,KTOP)
!1s   PTDIF(JL,KTOP)=ZTYY(JL,KTOP)*Z1BET(JL)
!1s   PQDIF(JL,KTOP)=ZQYY(JL,KTOP)*Z1BET(JL)
!1s ENDDO
!1s 
!1s !*         1.4    ELIMINATION FOR MIDDLE LAYERS.
!1s 
!1s DO JK=KTOP+1,KLEV-1
!1s   DO JL=KIDIA,KFDIA
!1s     ZGAM(JL,JK)=ZCC(JL,JK-1)*Z1BET(JL)
!1s     Z1BET(JL)=1.0_JPRB/(ZBB(JL,JK)-ZAA(JL,JK)*ZGAM(JL,JK))
!1s     PTDIF(JL,JK)=(ZTYY(JL,JK)-ZAA(JL,JK)*PTDIF(JL,JK-1))*Z1BET(JL)
!1s     PQDIF(JL,JK)=(ZQYY(JL,JK)-ZAA(JL,JK)*PQDIF(JL,JK-1))*Z1BET(JL)
!1s   ENDDO
!1s ENDDO
!1s 
!1s !*         1.5    BOTTOM LAYER, LINEAR RELATION BETWEEN LOWEST
!1s !                 MODEL LEVEL S AND Q AND FLUXES.
!1s 
!1s DO JL=KIDIA,KFDIA
!1s   ZGAM(JL,KLEV)=ZCC(JL,KLEV-1)*Z1BET(JL)
!1s   Z1BET(JL)=1.0_JPRB/(ZBB(JL,KLEV)-ZAA(JL,KLEV)*ZGAM(JL,KLEV))
!1s   ZBSL(JL)=(ZTYY(JL,KLEV)-ZAA(JL,KLEV)*PTDIF(JL,KLEV-1))*RVDIFTS*Z1BET(JL)
!1s   ZBQL(JL)=(ZQYY(JL,KLEV)-ZAA(JL,KLEV)*PQDIF(JL,KLEV-1))*RVDIFTS*Z1BET(JL)
!1s   ZQDP=1.0_JPRB/(PAPHM1(JL,KLEV)-PAPHM1(JL,KLEV-1))
!1s   ZASL(JL)=-RG*PTMST*ZQDP*RVDIFTS*Z1BET(JL)
!1s   ZAQL(JL)=ZASL(JL)
!1s ENDDO

!1s !*         1.0 - 1.5  OFFLINE SETUP

JK=KLEV
DO JL=KIDIA,KFDIA
    ZBSL(JL)=PSLM1(JL,JK)+PTMST*PSLGE(JL,JK)
    ZBQL(JL)=PQTM1(JL,JK)+PTMST*PQTE(JL,JK)
    ZASL(JL)=0.0_JPRB
    ZAQL(JL)=0.0_JPRB
ENDDO



!*         1.6    PREPARE ARRAY'S FOR CALL TO SURFACE ENERGY
!                 BALANCE ROUTINE

ZTSRF(KIDIA:KFDIA,1)=PSST(KIDIA:KFDIA)
ZTSRF(KIDIA:KFDIA,2)=PTICE(KIDIA:KFDIA)
ZTSRF(KIDIA:KFDIA,3)=PTSM1M(KIDIA:KFDIA)
ZTSRF(KIDIA:KFDIA,4)=PTSM1M(KIDIA:KFDIA)
ZTSRF(KIDIA:KFDIA,5)=PTSNOW(KIDIA:KFDIA)
ZTSRF(KIDIA:KFDIA,6)=PTSM1M(KIDIA:KFDIA)
ZTSRF(KIDIA:KFDIA,7)=PTSNOW(KIDIA:KFDIA)
ZTSRF(KIDIA:KFDIA,8)=PTSM1M(KIDIA:KFDIA)
IF (LEFLAKE) THEN
  DO JL=KIDIA,KFDIA
    IF(PHLICE(JL) > 1.E-9_JPRB) THEN ! 1.E-9 or H_ICE_MIN_FLK present
      ZTSRF(JL,9)=PTLICE(JL)
    ELSE
      ZTSRF(JL,9)=PTLWML(JL)
    ENDIF
  ENDDO
ENDIF

IF (LEURBAN) THEN
 ZTSRF(KIDIA:KFDIA,10)=PTSM1M(KIDIA:KFDIA)
 DO JL=KIDIA,KFDIA
  DO JK=1,10
    IF (ZTSRF(JL,JK) .NE. PTSRF(JL,JK) ) THEN
      WRITE(*,*) 'vdfdifh:',JL,JK,ZTSRF(JL,JK),PTSRF(JL,JK),ZTSRF(JL,JK)-PTSRF(JL,JK)
    ENDIF
  ENDDO
 ENDDO
ELSE
 DO JL=KIDIA,KFDIA
  DO JK=1,9
    IF (ZTSRF(JL,JK) .NE. PTSRF(JL,JK) ) THEN
      WRITE(*,*) 'vdfdifh:',JL,JK,ZTSRF(JL,JK),PTSRF(JL,JK),ZTSRF(JL,JK)-PTSRF(JL,JK)
    ENDIF
  ENDDO
 ENDDO
ENDIF

ZCOEF1=1.0_JPRB/(RG*RVDIFTS*PTMST)
DO JT=1,KTILES
  DO JL=KIDIA,KFDIA
    ZRHOCHU(JL,JT)=PCFHTI(JL,JT)*ZCOEF1
    ZRHOCQU(JL,JT)=PCFQTI(JL,JT)*ZCOEF1
  ENDDO
ENDDO

!*         1.7    CALL TO SURFACE ENERGY BALANCE ROUTINE
!                 REMEMBER: OUTPUT IS EXTRAPOLATED IN TIME
!! OFFLINE LIM IS NOT AVAILABLE, BUT VARIABLES NEEDED FOR INTERFACE
LSICOUP=.FALSE.
LNEMOLIMTHK=.FALSE.
ZSNTICE(:) = 0.0_JPRB
ZTHKICE(:) = 273._JPRB


!* Unstressed low vegetation evapotranspiration 
!* + open water pot. evap over tile 1/2 assuming water == air temperature 
!! Surface temperature for unstress evap (water)

ZTSRFU(KIDIA:KFDIA,:) = PTSRF(KIDIA:KFDIA,:)
! do jl=kidia,kfdia
!   do jt=1,9
!     if (ZTSRF(jl,jt) /= PTSRF(jl,jt) ) then
!       print*,jl,jt,ZTSRF(jl,jt) ,PTSRF(jl,jt)
!       stop -1
!     endif
!   enddo
! enddo
    
DO JL=KIDIA,KFDIA
  ZTSRFU(JL,1) = PTM1(JL,KLEV)
ENDDO
CALL SURFSEB   (YDSURF=YDSURF,KIDIA=KIDIA,KFDIA=KFDIA,KLON=KLON,KTILES=KTILES,KTVL=KTVL,KTVH=KTVH,&
 & PSSDP2=PSSDP2,PTMST=PTMST,PSSKM1M=PCPTSTIU,PTSKM1M=PTSKTI,PQSKM1M=PQSTI,&
 & PDQSDT=PDQSTI,PRHOCHU=ZRHOCHU,PRHOCQU=ZRHOCQU,&
 & PALPHAL=PCAIRTIU,PALPHAS=PCSATTIU,&
 & PSSRFL=PSSRFLTI,PFRTI=PFRTI,PTSRF=ZTSRFU,PLAMSK=PLAMSK,&
 & PSNM=PSNM,PRSN=PRSN,PHLICE=PHLICE, &
 & PSLRFL=PSLRFL,PTSKRAD=PTSKRAD,PEMIS=PEMIS,&
 & PASL=ZASL,PBSL=ZBSL,PAQL=ZAQL,PBQL=ZBQL,&
 & PTHKICE=ZTHKICE,PSNTICE=ZSNTICE,&
 !out
 & PJS=ZJS,PJQ=PJQU,PSSK=ZSSK,PTSK=ZTSK,&
 & PSSH=PSSH,PSLH=PSLH,PSTR=PSTR,PG0=PG0,&
 & PSL=ZSL,PQL=ZQL,&
 & LNEMOLIMTHK=LNEMOLIMTHK)  


!* fully Unstressed evapotranspiration 
DO JT=1,KTILES
  DO JL=KIDIA,KFDIA
    ZPCAIRTIU(JL,JT) = 1.0_JPRB
    ZPCSATTIU(JL,JT) = 1.0_JPRB
    ZPCPTSTIU(JL,JT) = PTSKTI(JL,JT)*RCPD*(1.0_JPRB+RVTMP2*&
   & (ZPCSATTIU(JL,JT)*PQSTI(JL,JT)+(1.0_JPRB-ZPCAIRTIU(JL,JT))*PQM1(JL,KLEV)))
  ENDDO
ENDDO

CALL SURFSEB   (YDSURF=YDSURF,KIDIA=KIDIA,KFDIA=KFDIA,KLON=KLON,KTILES=KTILES,KTVL=KTVL,KTVH=KTVH,&
 & PSSDP2=PSSDP2,PTMST=PTMST,PSSKM1M=ZPCPTSTIU,PTSKM1M=PTSKTI,PQSKM1M=PQSTI,&
 & PDQSDT=PDQSTI,PRHOCHU=ZRHOCHU,PRHOCQU=ZRHOCQU,&
 & PALPHAL=ZPCAIRTIU,PALPHAS=ZPCSATTIU,&
 & PSSRFL=PSSRFLTI,PFRTI=PFRTI,PTSRF=PTSRF,PLAMSK=PLAMSK,&
 & PSNM=PSNM,PRSN=PRSN,PHLICE=PHLICE, &
 & PSLRFL=PSLRFL,PTSKRAD=PTSKRAD,PEMIS=PEMIS,&
 & PASL=ZASL,PBSL=ZBSL,PAQL=ZAQL,PBQL=ZBQL,&
 & PTHKICE=ZTHKICE,PSNTICE=ZSNTICE,&
 !out
 & PJS=ZJS,PJQ=ZPJQU,PSSK=ZSSK,PTSK=ZTSK,&
 & PSSH=PSSH,PSLH=PSLH,PSTR=PSTR,PG0=PG0,&
 & PSL=ZSL,PQL=ZQL,&
 & LNEMOLIMTHK=LNEMOLIMTHK)  

! Swap memory to PJQU
DO JL=KIDIA,KFDIA
 PJQU(JL,2) = 0
 DO JT=1,KTILES
    PJQU(JL,2) = PJQU(JL,2)+ PFRTI(JL,JT)*ZPJQU(JL,JT)
 ENDDO
ENDDO


!* Normal call to surface energy balance ! 

CALL SURFSEB   (YDSURF=YDSURF,KIDIA=KIDIA,KFDIA=KFDIA,KLON=KLON,KTILES=KTILES,KTVL=KTVL,KTVH=KTVH,&
 & PSSDP2=PSSDP2,PTMST=PTMST,PSSKM1M=PCPTSTI,PTSKM1M=PTSKTI,PQSKM1M=PQSTI,&
 & PDQSDT=PDQSTI,PRHOCHU=ZRHOCHU,PRHOCQU=ZRHOCQU,&
 & PALPHAL=PCAIRTI,PALPHAS=PCSATTI,&
 & PSSRFL=PSSRFLTI,PFRTI=PFRTI,PTSRF=PTSRF,PLAMSK=PLAMSK,&
 & PSNM=PSNM,PRSN=PRSN,PHLICE=PHLICE, &
 & PSLRFL=PSLRFL,PTSKRAD=PTSKRAD,PEMIS=PEMIS,&
 & PASL=ZASL,PBSL=ZBSL,PAQL=ZAQL,PBQL=ZBQL,&
 & PTHKICE=ZTHKICE,PSNTICE=ZSNTICE,&
 !out
 & PJS=ZJS,PJQ=PJQ,PSSK=ZSSK,PTSK=ZTSK,&
 & PSSH=PSSH,PSLH=PSLH,PSTR=PSTR,PG0=PG0,&
 & PSL=ZSL,PQL=ZQL,&
 & LNEMOLIMTHK=LNEMOLIMTHK)
 
! Change PJQU (unstressed evap of tile 9 to match the actual evap of tile 9)
! this can be used as proxy for potential open water evaporation in "lake world" 
! to extract water from floodplains in CaMa-floodplains
! We cannot use direclty PJQ (PEVAPTI) because it get's multiplied by tile fraction latter
DO JL=KIDIA,KFDIA
  PJQU(JL,9) = PJQ(JL,9)
ENDDO

! 2nd call to surface energy balance 
! in case of skin temperature above freezing point 
! only active with the new multi-layer snow 
IF (LESNML) THEN
  ZTSRF(KIDIA:KFDIA,:) = PTSRF(KIDIA:KFDIA,:)
  ZLAMSK(KIDIA:KFDIA,:) = PLAMSK(KIDIA:KFDIA,:)
  DO JL=KIDIA,KFDIA
    IF (PFRTI(JL,5) >= 1.E-6_JPRB ) THEN
      IF (ZTSK(JL,5) >=RTT-1.E-6_JPRB ) THEN
        ZTSRF(JL,5)=RTT
        IF (PCIL(JL) >= 1._JPRB - 1.E-6_JPRB) THEN
          ZLAMSK(JL,5)= 1.E6_JPRB
        ELSE
          ZLAMSK(JL,5)= 100._JPRB
        ENDIF
      ENDIF
    ENDIF
    IF (PFRTI(JL,2) >= 1.E-6_JPRB ) THEN
      IF (ZTSK(JL,2) >= RTT-1.E-6_JPRB ) THEN
        ZTSRF(JL,2)=RTT
        ZLAMSK(JL,2)=1.E6_JPRB
      ENDIF
    ENDIF
  ENDDO
  ! call again surfseb replacing ZTSRF and ZLAMSK to avoid overshoot of tile 5 and tile 2
  ! skin temperatures
  CALL SURFSEB   (YDSURF=YDSURF,KIDIA=KIDIA,KFDIA=KFDIA,KLON=KLON,KTILES=KTILES,KTVL=KTVL,KTVH=KTVH,&
  & PSSDP2=PSSDP2, PTMST=PTMST,PSSKM1M=PCPTSTI,PTSKM1M=PTSKTI,PQSKM1M=PQSTI,&
  & PDQSDT=PDQSTI,PRHOCHU=ZRHOCHU,PRHOCQU=ZRHOCQU,&
  & PALPHAL=PCAIRTI,PALPHAS=PCSATTI,&
  & PSSRFL=PSSRFLTI,PFRTI=PFRTI,PTSRF=ZTSRF,PLAMSK=ZLAMSK,&
  & PSNM=PSNM,PRSN=PRSN,PHLICE=PHLICE, &
  & PSLRFL=PSLRFL,PTSKRAD=PTSKRAD,PEMIS=PEMIS,&
  & PASL=ZASL,PBSL=ZBSL,PAQL=ZAQL,PBQL=ZBQL,&
  & PTHKICE=ZTHKICE,PSNTICE=ZSNTICE,&
  !out
  & PJS=ZJS,PJQ=PJQ,PSSK=ZSSK,PTSK=ZTSK,&
  & PSSH=PSSH,PSLH=PSLH,PSTR=PSTR,PG0=PG0,&
  & PSL=ZSL,PQL=ZQL,&
  & LNEMOLIMTHK=LNEMOLIMTHK)
ENDIF 

! !ZCOEF2=0._JPRB
! do jt=1,KTILES
!  write(*,'(''XXjtile=''i4,1X,11(1X,f18.9))')jt,PFRTI(1,jt),PSSRFLTI(1,jt),PSLRFL(1),ZTSRF(1,jt),PTSKTI(1,jt),ZTSK(1,jt),pssh(1,jt),PSLH(1,jt),pstr(1,jt),PG0(1,jt),ZLAMSK(1,jt)
! ! ZCOEF2=ZCOEF2+PFRTI(1,jt)*pstr(1,jt)
! enddo

!*         1.7a   ADD SNOW EVAPORATION TO FLUXES

PJQ (KIDIA:KFDIA,7)=PJQ (KIDIA:KFDIA,7)+PEVAPSNW(KIDIA:KFDIA)
PSLH(KIDIA:KFDIA,7)=PSLH(KIDIA:KFDIA,7)+PEVAPSNW(KIDIA:KFDIA)*RLSTT

!*         1.7b   Flux boundary condition for 1D model (fluxes in W/m2)
!                 (Over-write output of SURFSEB)

IF (LDFLUX1D) THEN
  DO JT=1,KTILES
    DO JL=KIDIA,KFDIA
      ZJS(JL,JT)=PFSH1D(JL)+RCPD*PTSKTI(JL,JT)*RVTMP2*PFLH1D(JL)/RLVTT
      PJQ(JL,JT)=PFLH1D(JL)/RLVTT

      ZSSK(JL,JT)=ZBSL(JL)-ZJS(JL,JT)*(ZASL(JL)-1.0_JPRB/ZRHOCHU(JL,JT)) 
      ZTSK(JL,JT)=ZSSK(JL,JT)/(RCPD*(1._JPRB+RVTMP2*PQSTI(JL,JT)))
      PSSH(JL,JT)=PFSH1D(JL)
      PSLH(JL,JT)=PFLH1D(JL)
      PSTR(JL,JT)=PSLRFL(JL)
      PG0 (JL,JT)=PFSH1D(JL)+PFLH1D(JL)+PSLRFL(JL)+PSSRFLTI(JL,JT)
    ENDDO
  ENDDO
  DO JL=KIDIA,KFDIA
    ZSL(JL)=ZJS(JL,1)*ZASL(JL)+ZBSL(JL)
    ZQL(JL)=PJQ(JL,1)*ZAQL(JL)+ZBQL(JL)
  ENDDO
ENDIF

!*         1.7c   COMPUTE PARAMETERS AT NEW TIME LEVEL 

DO JT=1,KTILES
  DO JL=KIDIA,KFDIA
    PTSKTIP1(JL,JT)=ZTPFAC2*ZTSK(JL,JT)+ZTPFAC3*PTSKTI(JL,JT)
    ZQSP1=PQSTI(JL,JT)+PDQSTI(JL,JT)*(ZTSK(JL,JT)-PTSKTI(JL,JT))
  ENDDO
ENDDO

!*         1.7d   COPY LOWEST MODEL SOLUTION FROM SURFSEB

ZCOEF2=1.0_JPRB/RVDIFTS
DO JL=KIDIA,KFDIA
  PTDIF(JL,KLEV)=ZSL(JL)*ZCOEF2
  PQDIF(JL,KLEV)=ZQL(JL)*ZCOEF2
ENDDO

! !*         1.8    BACK-SUBSTITUTION.
! 
! DO JK=KLEV-1,KTOP,-1
!   DO JL=KIDIA,KFDIA
!     PTDIF(JL,JK)=PTDIF(JL,JK)-ZGAM(JL,JK+1)*PTDIF(JL,JK+1)
!     PQDIF(JL,JK)=PQDIF(JL,JK)-ZGAM(JL,JK+1)*PQDIF(JL,JK+1)
!   ENDDO
! ENDDO

  
IF (LHOOK) CALL DR_HOOK('VDFDIFH1S',1,ZHOOK_HANDLE)
END SUBROUTINE VDFDIFH1S
